#!/bin/csh -f
#
# aqpc: Shell script to run Aqua (complete analysis for all restraint types)
# in conjunction with, or separate from, Procheck-NMR.
# ===========================================================
# Ton Rullmann, Aug 18, 1997 - Aqua 1.4
# ===========================================================
# The script is set up analogous to procheck.scr.
# The actual work is done by calling the 'aqdrst' script.
# The restraint files should be in the current directory.
#
# This script expects at least one argument: the coordinates file.
#
# Optionally various arguments can be given, cf. the usage statement.
#
# The variable 'aquaroot' must exist and point to the Aqua top directory.
# The 'aqsetup' script should reside in that directory.
#

if ( $#argv < 1 ) goto usage
goto begin

#####################
# usage block
#####################
usage:
    echo 'usage: aqpc  [options]  coordinates_file  [ranges_files]'
    echo ''
    echo '   options (can be in any order):'
    echo '   [-c]  [-nosum]  [-r name]  [-ext]  [-multi] [-terseness value]'
    echo '   [-xplor|molmol]  [-r6ave|r6sum value]'
    echo ''
    echo '   -c:            prevent cleanup'
    echo '   -nosum:        do not make summary file'
    echo '   -r:            file name (minus extension) of Aqua restraint files'
    echo '                  (if not given, name is derived from coordinates file name)'
    echo '   -ext:          use the external (original) names to match restraint atoms'
    echo '                  to coordinate atoms (default is to use internal names)'
    echo '   -multi:        do not delete/combine similar restraints'
    echo '   -terseness value: [0-9] integer specifying decreasing verbosity'
    echo '   -xplor:        coordinates file is an XPLOR-style PDB file'
    echo '   -molmol:       coordinates file is a MOLMOL-style PDB file'
    echo '   -r6ave:        use R**(-6) averaging for all NOE restraints'
    echo '   -r6sum value:  use R**(-6) summing for all NOE restraints'
    echo '                  value is the number of monomers'
    exit (1)

#####################
# initialization
#####################
begin:
    set errstatus   = 0
    set cleanup     = 1
    set coordinp    = 0
    set restrname   = ""
    set stereo      = ""
    set nmono       = ""
    set oldmap      = "off"
    set multi       = 1
    set dosum       = 1
    set ident_all   = "pc"
    set ident_ste   = "pcs"
    set ident_rng   = "rng"
    set ident_rngs  = "rngs"
    set script_log  = "aqpc.log"
    setenv AQUATERSENESS 0     # most verbose messages

#####################
# process the options
#####################
options:
switch ( "$1" )
    case "-c":
	set cleanup = 0
	shift
	breaksw
    case "-r":
	if ( $#argv < 1 ) goto usage
	shift
	set restrname = $1
	shift
	breaksw
    case "-ext":
	set oldmap = "on"
	shift
	breaksw
    case "-multi":
	set multi = 0
	shift
	breaksw
    case "-nosum":
	set dosum = 0
	shift
	breaksw
    case "-xplor":
	set coordinp = 4
	shift
	breaksw
    case "-molmol":
	set coordinp = 5
	shift
	breaksw
    case "-r6ave":
	set stereo = "average"
	shift
	breaksw
    case "-r6sum":
	set stereo = "sum"
	shift
	set nmono = $1
	shift
	breaksw
    case "-terseness":
	shift        
	setenv AQUATERSENESS $1 # Inherited by all processes because using setenv
	shift
	breaksw
    default:
	goto endoptions
endsw
goto options

endoptions:
    if ( ( $#argv < 1 ) || ( $#argv > 2 ) ) goto usage

# Check options for acceptable ranges
if ( ( $AQUATERSENESS < 0 ) || ( $AQUATERSENESS >9 ) ) then
    echo "The parameter terseness (specified to be $AQUATERSENESS) can only"
    echo "take integer values from 0 through 9"
    exit(1)
endif

# Set level of terseness as a category which can be modified later on
if ( $AQUATERSENESS < 3 ) then
    set terselevel = "verbose"
else if ( $AQUATERSENESS < 6 ) then
    set terselevel = "curbed"
else
    set terselevel = "terse"
endif

# get started

if ( $terselevel == "verbose" ) then
    echo 'Aqua program directory:' $aquaroot
endif
source $aquaroot/aqsetup

if ( $terselevel == "verbose" ) then
    echo 'AQPC - Processing distance restraints'
endif

# set coordinates file
    set coords  = $1
    set coordsf = $1
    set coordsc = $1
    if ( `echo $coords | cut -c1 - ` == "%" ) then
	set coordsf = `echo $coords | cut -c2- - `
	set coordsc = `gawk 'NF>0{print $1;exit}' $coordsf`
    endif
    if ( ! -e $coordsf ) then
        echo '** AQPC - File' $coordsf 'does not exist'
        exit (1)
    endif
# set ranges file and idents
    set ranges = ""
    set ranges2 = ""
    if ( $#argv > 1 ) then
	set ranges = $2
	if ( ! -e $ranges ) then
	    echo '** AQPC - File' $ranges 'does not exist'
	    exit (1)
	endif
# add the ranges_ident; set the output ident
	set ranges = ( $ranges $ident_rng )
	set ident = $ident_rng
    else
	set ident = $ident_all
    endif
# identN is output ident for NOE only
    set identN = $ident
# stereo option; (re)set the output ident and the 'aqdrst make' command
    if ( $stereo == "" ) then
	set make = "make"
    else
	set make = "make_nodist"
	set identN = $ident_ste
	if ( "$ranges" != "" ) then
	    set ranges2 = ( $ranges $ident_rngs )
	    set identN = $ident_rngs
	endif
    endif
# setup project name
    set proj = `qgetproj $coordsf`
    if ( $terselevel == "verbose" ) then
        echo 'AQPC - Project name is' $proj
    endif
    qproject $proj
# set name (minus extension) of Aqua restraint file
    if ( $restrname == "" ) set restrname = $proj
# do we have to do anything?
    if ( ! -e $restrname.noe && ! -e $restrname.hb && ! -e $restrname.ss && ! -e $restrname.dis ) then
	echo '** AQPC - No AQUA distance restraint files ' $restrname'.? found'
        exit (1)
    endif
# remove old output files
    if ( -e $proj.nrv ) \rm $proj.nrv
    if ( -e $proj.nrc ) \rm $proj.nrc
    if ( -e $proj.nsm ) \rm $proj.nsm
    if ( -e $proj.hrv ) \rm $proj.hrv
    if ( -e $proj.hrc ) \rm $proj.hrc
    if ( -e $proj.hsm ) \rm $proj.hsm
    if ( -e $proj.srv ) \rm $proj.srv
    if ( -e $proj.src ) \rm $proj.src
    if ( -e $proj.ssm ) \rm $proj.ssm
    if ( -e $proj.drv ) \rm $proj.drv
    if ( -e $proj.drc ) \rm $proj.drc
    if ( -e $proj.dsm ) \rm $proj.dsm
# initialize the log file    
    \rm $script_log >& /dev/null
    
    if ( $terselevel != "verbose" ) then
        echo 'USE A LOWER VALUE FOR TERSENESS TO GET MORE OUTPUT IN THIS FILE IF NEEDED'    >> $script_log
        echo ' '                                                                            >> $script_log
    else
        echo 'AQPC - Processing distance restraints '           >> $script_log
        echo '--------------------------------------'           >> $script_log
        echo ' '                                                >> $script_log
        echo 'AQPC - Aqua program directory:' $aquaroot         >> $script_log
        echo ' '                                                >> $script_log
    endif
# make or update the aqua.cmm file
# first read the required characteristics from the restraint file
    foreach x ( $restrname.noe $restrname.hb $restrname.ss $restrname.dis )
	if ( ! -e $x ) continue
	set fmtline = (`sed -n '2p;2q' $x`)
	set fmt = $fmtline[4]
	set mr = "off"
	if ( $#fmtline == 5 ) then
	    if ( $fmtline[5] == "mr" ) set mr = "on"
	endif
	break
    end
# determine the coordinate file type
    if ( ! $coordinp ) then
	qguessc $coordsc >! qguessc.log
	set sts = $status
        if ( $terselevel == "verbose" ) then
	    cat qguessc.log
	    echo '*'        >> $script_log            
	    cat qguessc.log >> $script_log
	    echo '*'        >> $script_log            
        endif
	\rm qguessc.log
	if ( $sts == 0 ) then
	    echo "* ERROR: could not determine coordinate file type"
	    exit (1)
	else
	    set coordinp = $sts
	endif
    else
	if ( $coordinp == 4 ) then
	    echo "* Coordinate file type forced to PDB-XPLOR"
	    echo "* Coordinate file type forced to PDB-XPLOR" >> $script_log
	else if ( $coordinp == 5 ) then
	    echo "* Coordinate file type forced to PDB-MOLMOL"
	    echo "* Coordinate file type forced to PDB-MOLMOL" >> $script_log
	endif
    endif
# make or modify the aqua.cmm file
    touch aqua.cmm
    gawk -f $AQUASCRIPTSDIR/makecmm \
        -v coordinp=$coordinp -v restrinp=$fmt \
	-v oldmap=$oldmap     -v mr=0 \
        -v multi=$multi       -v terseness=$AQUATERSENESS \
        aqua.cmm > aqua.cmm.tmp
    \mv aqua.cmm.tmp aqua.cmm
# NOE restraints -- only these are done with stereo mode, if specified
    set type = "noe"
    set restrs = $restrname.$type
    if ( ! -e $restrs ) then
#	echo 'AQPC - File' $restrs 'does not exist'
    else
        if ( $terselevel == "verbose" ) then
	    echo ' '
            echo 'AQPC - Processing NOE restraints from' $restrs
        endif
	echo '********'                                 >> $script_log
        echo '** AQPC - Processing NOE restraints **'   >> $script_log
        echo '********'                                 >> $script_log
        aqdrst $make $ident_all $coords $restrs $type $ranges > aqdrst.log
	set errstatus = $status
	grep '^ *\*' aqdrst.log 
	cat aqdrst.log >> $script_log
	\rm aqdrst.log
	if ( ! $errstatus && $stereo != "" ) then
	    aqdrst stereo $stereo $nmono $ident_all noe $ident_ste $ranges2 > aqdrst.log
	    set errstatus2 = $status
	    grep '^ *\*' aqdrst.log 
	    cat aqdrst.log >> $script_log
	    \rm aqdrst.log
	    if ( $errstatus2 == 1 ) goto endaqpc
	endif
	echo ' ' >> $script_log
        if ( $terselevel == "verbose" ) then
	    echo "TIME0 aqdrst" >> $script_log
	    time                >> $script_log
	    echo ' '            >> $script_log
        endif
	if ( -e $proj'_'$type'_'$identN'.rvio' ) \mv $proj'_'$type'_'$identN'.rvio' $proj.nrv
	if ( -e $proj'_'$type'_'$identN'.rcnt' ) \mv $proj'_'$type'_'$identN'.rcnt' $proj.nrc
        if ( $terselevel == "verbose" ) then
            echo ' '                                                                >> $script_log
            echo 'AQPC - Output files have been renamed to:' $proj.nrv $proj.nrc    >> $script_log
        endif
	if ( $errstatus == 1 ) goto endaqpc
	if ( $errstatus != 10 && $dosum ) then
	    qsumm $proj.nrv $proj.nrc > $proj.nsm
	    echo ' '                                           >> $script_log
            if ( $terselevel == "verbose" ) then
	        echo "TIME0 qsumm "                            >> $script_log
	        time                                           >> $script_log
	        echo ' '                                       >> $script_log
	        echo 'AQPC - Summary of results on:' $proj.nsm >> $script_log
            endif
	endif
	set errstatus = 0
    endif
# H-BOND restraints
    set type = "hb"
    set restrs = $restrname.$type
    if ( ! -e $restrs ) then
#	echo 'AQPC - File' $restrs 'does not exist'
    else
        if ( $terselevel == "verbose" ) then
	    echo ' '
	    echo 'AQPC - Processing H-BOND restraints from' $restrs
        endif
	echo '*'                                                >> $script_log
        echo '********'                                         >> $script_log
        echo '** AQPC - Processing H-BOND restraints **'        >> $script_log
        echo '********'                                         >> $script_log
        aqdrst make $ident_all $coords $restrs $type $ranges >! aqdrst.log
	set errstatus = $status
	grep '^ *\*' aqdrst.log 
	cat aqdrst.log >> $script_log
	\rm aqdrst.log
	if ( -e $proj'_'$type'_'$ident'.rvio' ) \mv $proj'_'$type'_'$ident'.rvio' $proj.hrv
	if ( -e $proj'_'$type'_'$ident'.rcnt' ) \mv $proj'_'$type'_'$ident'.rcnt' $proj.hrc
        if ( $terselevel == "verbose" ) then
            echo ' '                                                            >> $script_log
	    echo 'AQPC - Output files have been renamed to:' $proj.hrv $proj.hrc    >> $script_log
        endif
	if ( $errstatus == 1 ) goto endaqpc
	if ( $errstatus != 10 && $dosum ) then
	    qsumm $proj.hrv $proj.hrc > $proj.hsm
            if ( $terselevel == "verbose" ) then
	        echo 'AQPC - Summary of results on:' $proj.hsm >> $script_log
            endif
	endif
	set errstatus = 0
    endif
# DISULPHIDE restraints
    set type = "ss"
    set restrs = $restrname.$type
    if ( ! -e $restrs ) then
#	echo 'AQPC - File' $restrs 'does not exist'
    else
        if ( $terselevel == "verbose" ) then    
            echo ' '
	    echo 'AQPC - Processing DISULPHIDE restraints from' $restrs
        endif
	echo '*'                                                >> $script_log
        echo '********'                                         >> $script_log
        echo '** AQPC - Processing DISULPHIDE restraints **'    >> $script_log
        echo '********'                                         >> $script_log
        aqdrst make $ident_all $coords $restrs $type $ranges >! aqdrst.log
	set errstatus = $status
	grep '^ *\*' aqdrst.log 
	cat aqdrst.log >> $script_log
	\rm aqdrst.log
	if ( -e $proj'_'$type'_'$ident'.rvio' ) \mv $proj'_'$type'_'$ident'.rvio' $proj.srv
	if ( -e $proj'_'$type'_'$ident'.rcnt' ) \mv $proj'_'$type'_'$ident'.rcnt' $proj.src
        if ( $terselevel == "verbose" ) then
	    echo ' '                                                             >> $script_log
            echo 'AQPC - Output files have been renamed to:' $proj.srv $proj.src >> $script_log
        endif
	if ( $errstatus == 1 ) goto endaqpc
	if ( $errstatus != 10 && $dosum ) then
	    qsumm $proj.srv $proj.src > $proj.ssm
            if ( $terselevel == "verbose" ) then
	        echo 'AQPC - Summary of results on:' $proj.ssm >> $script_log
            endif
	endif
	set errstatus = 0
    endif
# GENERIC DISTANCE restraints
    set type = "dis"
    set restrs = $restrname.$type
    if ( ! -e $restrs ) then
#	echo 'AQPC - File' $restrs 'does not exist'
    else
        if ( $terselevel == "verbose" ) then
	    echo ' '
	    echo 'AQPC - Processing GENERIC DISTANCE restraints from' $restrs
        endif
	echo '*'                                                    >> $script_log
        echo '********'                                             >> $script_log
        echo '** AQPC - Processing GENERIC DISTANCE restraints **'  >> $script_log
        echo '********'                                             >> $script_log
        aqdrst make $ident_all $coords $restrs $type $ranges >! aqdrst.log
	set errstatus = $status
	grep '^ *\*' aqdrst.log 
	cat aqdrst.log >> $script_log
	\rm aqdrst.log
	if ( -e $proj'_'$type'_'$ident'.rvio' ) \mv $proj'_'$type'_'$ident'.rvio' $proj.drv
	if ( -e $proj'_'$type'_'$ident'.rcnt' ) \mv $proj'_'$type'_'$ident'.rcnt' $proj.drc
        if ( $terselevel == "verbose" ) then
	    echo ' '                                                             >> $script_log
            echo 'AQPC - Output files have been renamed to:' $proj.drv $proj.drc >> $script_log
        endif
	if ( $errstatus == 1 ) goto endaqpc
	if ( $errstatus != 10 && $dosum ) then
	    qsumm $proj.drv $proj.drc > $proj.dsm
            if ( $terselevel == "verbose" ) then
	        echo 'AQPC - Summary of results on:' $proj.dsm >> $script_log
            endif
	endif
	set errstatus = 0
    endif
# TORSION ANGLE restraints
    set type = "tor"
    set restrs = $restrname.$type
    if ( ! -e $restrs ) then
#	echo 'AQPC - File' $restrs 'does not exist'
    else
	if ( "$restrs" != "$proj.$type" ) then
	    echo 'AQPC - Copying TORSION ANGLE restraints from' $restrs 'to' $proj.$type
	    \cp $restrs $proj.$type
	endif
    endif

# cleanup files, if not turned off
endaqpc:
    if ( $cleanup ) then
	qproject $proj
        if ( $terselevel == "verbose" ) then
            echo ' '
    	    echo 'AQPC - Cleaning up for project' $proj
        endif
	aqdrst clean $ident_all
	if ( $#argv > 1 ) then
	    aqdrst clean $ident_rng
	endif
	if ( $stereo != "" ) then
	    aqdrst clean $ident_ste
	    if ( $#argv > 1 ) then
		aqdrst clean $ident_rngs
	    endif
	endif
    endif
    
    if ( $dosum ) then
	set list = ( `ls $proj.*sm` )
    else
	set list = ""
    endif
    echo '********'                                         >> $script_log
    echo '** AQPC - Summary of results on:' $list '**'      >> $script_log
    echo '********'                                         >> $script_log
    if ( $terselevel == "verbose" ) then
        echo "AQPC - Log file: $script_log"
        echo 'AQPC - Summary of results on:' $list
    endif
    if ( $errstatus ) exit (1)
